function [rules_quad, rules_comb] = summarize_plots_for_slides(covariate,wave,cost,winter_prevuse,SMC,fbase)
% This function plots the treatment rules learned using the EWM method for
% slides, without marginal distributions

% Inputs: 
% (1) covariate: income, size, vintage
% (2) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0)
% (3) cost: =true if welfare measure is in terms of cost savings using private marginal cost 
% =false if welfare measure is in terms of kWh reduction
% (4) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (5) SMC: =1 use social marginal cost; =0 use retail electricity price
% (6) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) rules_quad: only quadrant rule
% (2) rules_comb: quadrant rule overlaid with cubic rule

%% input treatment rule results
% input quadrant results
if (cost)
    tcost = 0.765; 
else
    tcost = 0;
end

switch winter_prevuse 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end
filename_coefs=sprintf('coef_quadrant_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
file = sprintf('%s_%s',fbase,filename_coefs);
S_quadrant=load(file);

%extract content of S_quadrant
c_fieldnames = fieldnames(S_quadrant);
for ifield = 1:length(c_fieldnames)
    field_ = c_fieldnames{ifield};
    eval(sprintf('%s=S_quadrant.%s;',field_,field_))
end

% input cubic results
switch winter_prevuse 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end
filename_coefs=sprintf('coef_cubic_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
file = sprintf('%s_%s',fbase,filename_coefs);
S=load(file);
%extract content of S
c_fieldnames = fieldnames(S);
for ifield = 1:length(c_fieldnames)
    field_ = c_fieldnames{ifield};
    eval(sprintf('%s=S.%s;',field_,field_))
end

% "beta" conflicts with a built-in function, have to do this extra
% assignment
beta=S.beta;

%% Generate plot inputs 
if min_sign1<0 && min_sign2<0
    x_coord=[min(v_x1) min(v_x1) v_x1(min_i1) v_x1(min_i1)];
    y_coord=[min(v_x2) v_x2(min_i2) v_x2(min_i2) min(v_x2)];
elseif min_sign1<0 && min_sign2>0 
    x_coord=[min(v_x1) min(v_x1) v_x1(min_i1) v_x1(min_i1)];
    y_coord=[v_x2(min_i2) max(v_x2) max(v_x2) v_x2(min_i2)];
elseif min_sign1>0 && min_sign2<0 
    x_coord=[v_x1(min_i1) v_x1(min_i1) max(v_x1) max(v_x1)];
    y_coord=[min(v_x2) v_x2(min_i2) v_x2(min_i2) min(v_x2)];
elseif min_sign1>0 && min_sign2>0    
    x_coord=[v_x1(min_i1) v_x1(min_i1) max(v_x1) max(v_x1)];
    y_coord=[v_x2(min_i2) max(v_x2) max(v_x2) v_x2(min_i2)];
end 

% cubic plot inputs 
line_x1=[min(x1_wave):1:max(x1_wave)]';
% calculate the usage cutoff for treatment at different income levels 
line_usage=Xscale(1,4)*...
               (-(beta(1)+line_x1*beta(2)./Xscale(1,1)...
               +(line_x1.^2)*beta(3)./Xscale(1,2)...
               +(line_x1.^3)*beta(4)./Xscale(1,3))./beta(5));
           
if (beta(5)>0) % treatment rule is increasing in pre-treatment consumption
    select=(line_usage<=max(prevuse_wave));
    patch_x=[line_x1(select); flipud(line_x1(select))];
    patch_y=[max(line_usage(select),min(prevuse_wave)); max(prevuse_wave)*ones(sum(select),1) ];
else           % treatment rule is decreasing in pre-treatment consumption
    select=(line_usage>=min(prevuse_wave));
    patch_x=[line_x1(select); flipud(line_x1(select))];
    patch_y=[min(line_usage(select),max(prevuse_wave)); min(prevuse_wave)*ones(sum(select),1) ];
end

%% Quadrant plot only 
x0=10;
y0=10;
width=800;
height=600;

rules_quad=figure; hold on;
set(gcf,'position',[x0,y0,width,height]);

patch (x_coord, y_coord, 'red', 'LineStyle', '-', 'FaceAlpha', 0.3,'EdgeColor', 'none')
if covariate=="vintage"
    axis([min(patch_x),max(patch_x),min(v_x2),max(v_x2)]);
else
    axis([min(v_x1),max(v_x1),min(v_x2),max(v_x2)]);
end
ylabel('Pretreatment Usage [kWh/mo]')
if covariate=="income"
    xlabel('Income [$ k]');
elseif covariate=="size"
    xlabel('House size [sq ft]');
elseif covariate=="vintage"
    xlabel('House construction year');
end
legend('Quadrant rule treat')
set(findall(gcf,'-property','FontSize'),'FontSize',24)

hold off;

%% Quadrant and cubic rules with no density subplots
x0=10;
y0=10;
width=800;
height=600;

rules_comb=figure; hold on;
set(gcf,'position',[x0,y0,width,height]);

patch (x_coord, y_coord, 'red', 'LineStyle', '-', 'FaceAlpha', 0.3,'EdgeColor', 'none')
patch (patch_x,patch_y,'blue', 'LineStyle', '-', 'FaceAlpha', 0.45,'EdgeColor', 'none')
if covariate=="vintage"
    axis([min(patch_x),max(patch_x),min(v_x2),max(v_x2)]);
else
    axis([min(v_x1),max(v_x1),min(v_x2),max(v_x2)]);
end
ylabel('Pretreatment Usage [kWh/mo]')
if covariate=="income"
    xlabel('Income [$ k]');
elseif covariate=="size"
    xlabel('House size [sq ft]');
elseif covariate=="vintage"
    xlabel('House construction year');
end
legend('Quadrant rule treat','Cubic rule treat')
set(findall(gcf,'-property','FontSize'),'FontSize',24)

end